/*  expupto.cpp  by K.Tsuru  */

/************************************
One of auxiliary tools for SN library
*************************************/
#ifndef  DEF_CONST_H
#include "defconst.h"
#endif  // DEF_CONST_H

#ifndef TOOLS_H
#include "tools.h"
#endif //TOOLS_H

// It returns the number of figures n!.
static double Stirling(double n){
    static const double C1 = 0.5 * log10(2 * M_PI) +1.0, C2 = log10(M_E);
    return C1 - C2 * n + (n + 0.5) * log10(n);
}

static double fig, log10X;

static double fx(double n) {
    return fig + n * log10X - Stirling(n);
}
/**************************************************
It returns how many terms needs to evaluate exp(x) in a
given precision.
 x^n
----- = 10^(-f)
 n!
i.e.
 f + n log10(x) - log10(n!) = 0
parameter log10x : Prease give log10(x).
***************************************************/

long upToExpSeries(long precision, double log10x) {
    fig = precision;    log10X = log10x;
    if ( fx(1.0) < 0) return 1L;    // x < 10^(-f)
    return (long)(bisection(1, precision, 1.0e-5, fx) + 1.0);
}
/// a test program
#if 0
static long upToExp1(long precision) {
    return upToExpSeries(precision, 0.0);
}

int main() {
    long L;
    for (L = 500; L < 5000; L += 500) {
        long p = upToExp1(L);
        printf("%ld : %g, p = %ld, %g, %g\n", L, Stirling((double)p), p, fx(p-1), fx(p));
    }
    L = 1000;
    for (double x = 1.0; x > 1.0e-200; x /= 1.0e10) {
        long p = upToExpSeries(L, log10(x));
        printf("%ld : p = %ld, %g, %g (x = %g)\n", L, p, fx(p-1), fx(p), x);
    }
    return 0;
}
#endif